*diff in diff analysis by voucher redesign
*last modified: 18 April 2025
*last modified by: Kathya Tapia-Schythe
*this do file implements the eventstudyinteract package to estimate individual 2x2 doubly robust did estimates and aggregate to obtain an att 

*-------------------------------------------------------------------------------
*--- Preamble
*-------------------------------------------------------------------------------

clear all
set more off
set rmsg on

set maxvar 120000
set emptycells drop

*ssc install csdid, replace
*ssc install drdid, replace
*net install csdid2, from("https://friosavila.github.io/stpackages")

*-------------------------------------------------------------------------------
*--- Directories and Log 
*-------------------------------------------------------------------------------

local data_dir ./data/cleaned
local raw_dir ./data/raw
local out_dir ./analysis/output
local graph_dir ./analysis/output/graphs
local tab_dir ./analysis/output/tables
local log_dir ./documentation/logs
local date: display %tdYY-NN-DD date(c(current_date), "DMY")
di "`date'"
capture log close

log using `log_dir'/voucher_sa_tip`date', replace

*-------------------------------------------------------------------------------
*--- Data and Frames 
*-------------------------------------------------------------------------------

*import cleaned eWIC implementation info into own frame
*eWIC implementation variable is at county fips/fiscal year level 
*partition data into state/treatment year groups

*import cleaned ewic implementation info into own frame
mkf ewic_fy
frame ewic_fy{
	use ./data/cleaned/ewic_rollout.dta, clear
	rename year fiscalyear
	drop pre* post* 
}

mkf income
frame income{
	use ./data/cleaned/countyyearmedianincome.dta, clear
	keep if year == 2005
	drop year
	destring ct_fips, replace
}

mkf countychars
frame countychars{
	use ./data/cleaned/countystats.dta, clear
	keep if year == 2005
	drop year
	destring ct_fips, replace
}

mkf demos
frame demos{
	use ./data/cleaned/demos.dta, clear
	keep if year == 2005
	drop year
	destring ct_fips, replace
}

cwf ewic_fy
frame put ct_fips ev_year state, into(timings)

cwf timings
duplicates drop
encode state, generate(state_code)

frlink m:1 ct_fips, frame(countychars) generate(countylink)
frget medincome sharepoor shareltwic sharelt2fpl sharesnap sharecash sharecashorsnap, from(countylink)
frlink m:1 ct_fips, frame(demos) generate(demoslink)
frget total_pop under5_pop hispanic_total_pop total_black_pop region division urbanrural share_black share_hispanic share_under5, from(demoslink)
gen log_pop = log(total_pop)

*import TIP store data into own frame
mkf tip_fy
frame tip_fy{
	local data_dir ./data/cleaned
	use `data_dir'/tip_auth_sq, clear
	*drop "Direct Distribution Center" and "Home Food Delivery Contractor" types
	gen f = (vendor_type1 == 3 & (vendor_type2 == 4 | missing(vendor_type2))) | (vendor_type1 == 4 & (vendor_type2 == 3 | missing(vendor_type2)))
	bys tip_id: egen flag = mode(f)
	bys tip_id: replace vendor_type1 = vendor_type1[_n+1] if flag == 0 & f == 1
	drop if flag == 1
	drop flag f tip_year_id
	*fill in gaps in state fips using modes
	generate str_fips = string(ct_fips)
	replace str_fips = "0" + str_fips if strlen(str_fips) == 4
	generate st_fips = substr(str_fips, 1, 2)
	destring st_fips, replace 
	
	*drop Vermont and mississippi because of different regiemes prior to eWIC
	*drop NV because ewic turn on then off then on again and we don't have reliable timings
	*drop missing st_fips
	drop if st_fips == 50 | st_fips == 32 | st_fips == 28 | missing(st_fips)

	*flag stores that specialize in WIC
	gen a50 = vendor_type1 == 1 | vendor_type1 == 7 | vendor_type2 == 1 | vendor_type2 == 7
	
	*link in low income/high poverty info
	frlink m:1 ct_fips, frame(income) generate(inc_link)
	frget lowinc highpov highwicelig, from(inc_link)

	*link in controls from timings frame
	frlink m:1 ct_fips, frame(timings) generate(timings_link)
	frget total_pop under5_pop hispanic_total_pop total_black_pop region division urbanrural share_black share_hispanic share_under5 medincome sharepoor shareltwic sharelt2fpl sharesnap sharecash sharecashorsnap, from(timings_link)

	*link ewic implementation info
	frlink m:1 ct_fips fiscalyear, frame(ewic_fy) generate(ewic_link)
	frget ev_year , from(ewic_link)
	* assume treated at earliest exposure
	bys tip_id: egen ey_min = min(ev_year)
	
	* assume treated at earliest exposure
	replace ev_year = ey_min if ev_year != ey_min & !missing(ey_min)
	bys tip_id (ev_year): assert ev_year[1] == ev_year[_N]

	*assert no dups at future level of xtset 
	duplicates report tip_id fiscalyear 
	assert r(N) == r(unique_value)
	xtset tip_id fiscalyear
	compress
}

cwf tip_fy

*-------------------------------------------------------------------------------
*--- Variables for Sun & Abraham (2021)
*-------------------------------------------------------------------------------

foreach c in never last {

 cap drop redesigntime redesigngroup timeToTreat L*_beforegroup F*_beforegroup L*_aftergroup F*_aftergroup L*_beforetime F*_beforetime L*_aftertime F*_aftertime
 
 cap drop ebt*
	
 *Control: never-treated units
 if ("`c'"=="never") {
  tab ev_year, m
  gen ev_year_never=ev_year
  replace ev_year_never=. if ev_year_never>2018
  gen never_year=(ev_year_never==.)
  tab never_year, m
  
  *Group: Before and After Voucher Redesign
  gen redesigngroup=(ev_year>2009)
  tab redesigngroup, m

  *Time: Before and After Voucher Redesign
  gen redesigntime=(fiscalyear>2009)
  tab redesigntime, m
  
  *Group: Generate treatment
  gen ebtgroup_before=(fiscalyear>=ev_year_never & redesigngroup==0)
  tab ebtgroup_before, m
  gen ebtgroup_after= (fiscalyear>=ev_year_never & redesigngroup==1)
  tab ebtgroup_after, m
  
  *Time: Generate treatment
  gen ebttime_before=(fiscalyear>=ev_year_never & redesigntime==0)
  tab ebttime_before, m
  gen ebttime_after= (fiscalyear>=ev_year_never & redesigntime==1)
  tab ebttime_after, m

  *Relative time
  gen timeToTreat=fiscalyear-ev_year_never
  tab timeToTreat, m

 *Before
  forvalues l = 0/4 {
   gen L`l'_beforegroup = (timeToTreat==`l' & redesigngroup==0)
  }
  gen L5_beforegroup = (timeToTreat>=5 & redesigngroup==0)

  forvalues f = 1/4 {
   gen F`f'_beforegroup = (timeToTreat==-`f' & redesigngroup==0)
  }
  gen F5_beforegroup = (timeToTreat<=-5 & redesigngroup==0)

  drop F1_beforegroup // Baseline t=-1 
  gen F1_beforegroup = 0 

 *Before
  forvalues l = 0/4 {
   gen L`l'_beforetime = (timeToTreat==`l' & redesigntime==0)
  }
  gen L5_beforetime = (timeToTreat>=5 & redesigntime==0)

  forvalues f = 1/4 {
   gen F`f'_beforetime = (timeToTreat==-`f' & redesigntime==0)
  }
  gen F5_beforetime = (timeToTreat<=-5 & redesigntime==0)

  drop F1_beforetime // Baseline t=-1 
  gen F1_beforetime = 0 

  *After
  forvalues l = 0/4 {
   gen L`l'_aftergroup = (timeToTreat==`l' & redesigngroup==1)
  }
  gen L5_aftergroup = (timeToTreat>=5 & redesigngroup==1)

  forvalues f = 1/4 {
   gen F`f'_aftergroup = (timeToTreat==-`f' & redesigngroup==1)
  }
  gen F5_aftergroup = (timeToTreat<=-5 & redesigngroup==1)

  drop F1_aftergroup // Baseline t=-1 
  gen F1_aftergroup = 0 

  *After
  forvalues l = 0/4 {
   gen L`l'_aftertime = (timeToTreat==`l' & redesigntime==1)
  }
  gen L5_aftertime = (timeToTreat>=5 & redesigntime==1)

  forvalues f = 1/4 {
   gen F`f'_aftertime = (timeToTreat==-`f' & redesigntime==1)
  }
  gen F5_aftertime = (timeToTreat<=-5 & redesigntime==1)

  drop F1_aftertime // Baseline t=-1 
  gen F1_aftertime = 0 
  
  local never "_never"
 }

*Control: last-treated unit
 if ("`c'"=="last") {
  sum ev_year
  gen last_year=(ev_year==r(max))
  tab last_year, m
  
  *Group: Before and After Voucher Redesign
  gen redesigngroup=(ev_year>2009)
  tab redesigngroup, m

  *Time: Before and After Voucher Redesign
  gen redesigntime=(fiscalyear>2009)
  tab redesigntime, m
  
  *Group: Generate treatment
  gen ebtgroup_before=(fiscalyear>=ev_year & redesigngroup==0)
  replace ebtgroup_before=. if ev_year>2018 & ev_year<2023 // Only use last-treated
  tab ebtgroup_before, m
  gen ebtgroup_after= (fiscalyear>=ev_year & redesigngroup==1)
  replace ebtgroup_after=. if ev_year>2018 & ev_year<2023 // Only use last-treated
  tab ebtgroup_after, m
  
  *Time: Generate treatment
  gen ebttime_before=(fiscalyear>=ev_year & redesigntime==0)
  replace ebttime_before=. if ev_year>2018 & ev_year<2023 // Only use last-treated
  tab ebttime_before, m
  gen ebttime_after= (fiscalyear>=ev_year & redesigntime==1)
  replace ebttime_after=. if ev_year>2018 & ev_year<2023 // Only use last-treated
  tab ebttime_after, m

  *Relative time
  gen timeToTreat=fiscalyear-ev_year
  tab timeToTreat, m

 *Before
  forvalues l = 0/4 {
   gen L`l'_beforegroup = (timeToTreat==`l' & redesigngroup==0)
  }
  gen L5_beforegroup = (timeToTreat>=5 & redesigngroup==0)

  forvalues f = 1/4 {
   gen F`f'_beforegroup = (timeToTreat==-`f' & redesigngroup==0)
  }
  gen F5_beforegroup = (timeToTreat<=-5 & redesigngroup==0)

  drop F1_beforegroup // Baseline t=-1 
  gen F1_beforegroup = 0 

 *Before
  forvalues l = 0/4 {
   gen L`l'_beforetime = (timeToTreat==`l' & redesigntime==0)
  }
  gen L5_beforetime = (timeToTreat>=5 & redesigntime==0)

  forvalues f = 1/4 {
   gen F`f'_beforetime = (timeToTreat==-`f' & redesigntime==0)
  }
  gen F5_beforetime = (timeToTreat<=-5 & redesigntime==0)

  drop F1_beforetime // Baseline t=-1 
  gen F1_beforetime = 0 

  *After
  forvalues l = 0/4 {
   gen L`l'_aftergroup = (timeToTreat==`l' & redesigngroup==1)
  }
  gen L5_aftergroup = (timeToTreat>=5 & redesigngroup==1)

  forvalues f = 1/4 {
   gen F`f'_aftergroup = (timeToTreat==-`f' & redesigngroup==1)
  }
  gen F5_aftergroup = (timeToTreat<=-5 & redesigngroup==1)

  drop F1_aftergroup // Baseline t=-1 
  gen F1_aftergroup = 0 

  *After
  forvalues l = 0/4 {
   gen L`l'_aftertime = (timeToTreat==`l' & redesigntime==1)
  }
  gen L5_aftertime = (timeToTreat>=5 & redesigntime==1)

  forvalues f = 1/4 {
   gen F`f'_aftertime = (timeToTreat==-`f' & redesigntime==1)
  }
  gen F5_aftertime = (timeToTreat<=-5 & redesigntime==1)

  drop F1_aftertime // Baseline t=-1 
  gen F1_aftertime = 0 
  
  local never
 }
 
*-------------------------------------------------------------------------------
*--- Estimates: 2x2 All, Chain, Independent for Redesign Based on Timing of Event Year
*-------------------------------------------------------------------------------

 cd `out_dir'

 *Samples
 cap drop all indep
 gen all=1
 gen indep=(chain==0 & chain!=.)

 lab var ebtgroup_before "WIC EBT $times$ Before Package Change"
 lab var ebtgroup_after  "WIC EBT $times$ After Package Change"

 foreach sample in all chain indep {

  #delimit ;
  eventstudyinteract auth ebtgroup_before ebtgroup_after if `sample'==1, 
				  vce(cluster st_fips)  
                                  cohort(ev_year`never') 
				  control_cohort(`c'_year)
				  absorb(tip_id fiscalyear);
  #delimit cr
  
  mat beta = e(b_iw)
  mat sigma = e(V_iw)
  erepost b=beta V=sigma, rename
  estimates save `out_dir'/sa_voucher_group_`c'_`sample'.ster, replace

 }

*-------------------------------------------------------------------------------
*--- Estimates: Event Studies All, Chain, Independent for Redesign Based on Timing of Event Year
*-------------------------------------------------------------------------------

cd `out_dir'

local times F5_beforegroup F4_beforegroup F3_beforegroup F2_beforegroup L0_beforegroup ///
            L1_beforegroup L2_beforegroup L3_beforegroup L4_beforegroup L5_beforegroup ///
	    F5_aftergroup  F4_aftergroup  F3_aftergroup  F2_aftergroup  L0_aftergroup ///
            L1_aftergroup  L2_aftergroup  L3_aftergroup  L4_aftergroup  L5_aftergroup 

foreach sample in all chain indep {
	
 preserve
 
 *Periods
 gen t=.
 local i=1
 foreach n of numlist -4/4 {
  replace t=`n' in `i'
  local ++i
 }
 	
 #delimit ;
 eventstudyinteract auth `times' if `sample'==1, vce(cluster st_fips)  
                                 cohort(ev_year`never') 
				 control_cohort(`c'_year)
				 absorb(tip_id fiscalyear);
 #delimit cr	
 
 local df = e(df_r)
   
 mat b  = e(b_iw)
 mat V  = e(V_iw)
     
 *Estimates
 mat blead_before = b[1,2..4]
 mat blag_before  = b[1,5..9]
 mat beta_before = blead_before,0,blag_before
 mat beta_before = beta_before'
 svmat beta_before, names(beta_before)
 rename beta_before1 beta_before
 
 mat blead_after  = b[1,12..14]
 mat blag_after   = b[1,15..19]
 mat beta_after = blead_after,0,blag_after
 mat beta_after = beta_after'
 svmat beta_after, names(beta_after)
 rename beta_after1 beta_after
  
 mata st_matrix("se",sqrt(diagonal(st_matrix("e(V_iw)"))))
 mat se = se'
  
 mat slead_before = se[1,2..4]
 mat slag_before  = se[1,5..9]
 mat se_before = slead_before,0,slag_before
 mat se_before = se_before'
 svmat se_before, names(se_before)
 rename se_before1 se_before
 
 mat slead_after  = se[1,12..14]
 mat slag_after   = se[1,15..19]
 mat se_after = slead_after,0,slag_after
 mat se_after = se_after'
 svmat se_after, names(se_after)
 rename se_after1 se_after
 
 *95% CI
 gen lb_before = beta_before-invttail(`df',.025)*se_before 
 gen ub_before = beta_before+invttail(`df',.025)*se_before 
 
 gen lb_after = beta_after-invttail(`df',.025)*se_after  
 gen ub_after = beta_after+invttail(`df',.025)*se_after 
 
 ereturn post b V
 
 *Test 
 lincom (L0_before+L1_before+L2_before+L3_before+L4_before)/5 - ///
        (L0_after+L1_after+L2_after+L3_after+L4_after)/5

 local pv = round(r(p), 0.001)
	   
 foreach r in before after {
 
  gen att_`r'=.
  gen latt_`r'=.
  gen uatt_`r'=.

 *ATT pre period
 lincom (F4_`r' + F3_`r' + F2_`r')/3
 replace att_`r'=r(estimate) if t<-1 & t!=.
 replace latt_`r'=r(lb) if t<-1 & t!=.
 replace uatt_`r'=r(ub) if t<-1 & t!=.

 *ATT post period
 lincom (L0_`r' + L1_`r' + L2_`r' + L3_`r' + L4_`r')/5
 replace att_`r'=r(estimate) if t>=0 & t!=.
 replace latt_`r'=r(lb) if t>=0 & t!=.
 replace uatt_`r'=r(ub) if t>=0 & t!=.
 
 }

*-------------------------------------------------------------------------------
*--- Export Results: Figure
*-------------------------------------------------------------------------------

 rename t t_before
 gen t_after=t_before+.15
 
 
 if ("`sample'"=="all")   local color1 black
 if ("`sample'"=="chain") local color1 emerald
 if ("`sample'"=="indep") local color1 orange
 
 if ("`sample'"=="all")   local color2 gs5
 if ("`sample'"=="chain") local color2 dkgreen
 if ("`sample'"=="indep") local color2 dkorange
  
 if ("`sample'"=="all")   local msym1 O
 if ("`sample'"=="chain") local msym1 O
 if ("`sample'"=="indep") local msym1 T
 
 if ("`sample'"=="all")   local msym2 Oh
 if ("`sample'"=="chain") local msym2 Oh
 if ("`sample'"=="indep") local msym2 Th
 
 sum att_before if t_after>=0
 local earlyatt = round(r(mean), 0.001)
 sum att_after if t_after>=0
 local lateatt = round(r(mean), 0.001)
 
 *ssc install coefplot, replace
 *ssc install blindschemes, replace 
 set scheme plotplainblind
 graph set window fontface "Times New Roman"

 #delimit ; 
 twoway rarea uatt_before latt_before t_before if t_before<0, 
	color(gs7%50)  ||
        rarea uatt_before latt_before t_before if t_before>=0, 
	color(gs7%50) || 
	scatter beta_before t_before, mcolor(`color1') msymbol(`msym1') 
        msize(medlarge) || 
	rcap ub_before lb_before t_before, msize(vtiny) lcolor(`color1') ||
	line att_before t_before if t_before<0, lcolor(`color1') 
	lpattern(solid) ||
	line att_before t_before if t_before>=0, lcolor(`color1') 
	lpattern(solid)	||
	rarea uatt_after latt_after t_after if t_after<0, 
	color(gs12%50)  ||
        rarea uatt_after latt_after t_after if t_after>=0, 
	color(gs12%50) || 
	scatter beta_after t_after, mcolor(`color2') msymbol(`msym2') 
        msize(medlarge) || 
	rcap ub_after lb_after t_after, msize(vtiny) lcolor(`color2') ||
	line att_after t_after if t_after<0, lcolor(`color2') 
	lpattern(solid) ||
	line att_after t_after if t_after>=0, lcolor(`color2') 
	lpattern(solid) 
	yline(0, lcolor(black) lpattern(dash)) 
	xline(0, lcolor(red) lpattern(dot) lwidth(medthick))
	ytitle("ATT estimate") xtitle("Event Year")
	note("Post-EBT ATT estimate for counties implementing before package change is  `earlyatt' and after" "package change is `lateatt'." "P-value for difference in their average effects is `pv'")
	legend(order(3 "Before Package Change" 
	             9 "After Package Change") rows(1) position(6));
 #delimit cr

 graph export "`graph_dir'/voucher_group_sa_tip_`c'_`sample'.png", replace

 restore
  
}

*-------------------------------------------------------------------------------
*--- Estimates: 2x2 All, Chain, Independent for Redesign Based on Timing of Calendar Year
*-------------------------------------------------------------------------------

 cd `out_dir'

 lab var ebttime_before "WIC EBT $times$ Before Package Change"
 lab var ebttime_after  "WIC EBT $times$ After Package Change"

 foreach sample in all chain indep {

  #delimit ;
  eventstudyinteract auth ebttime_before ebttime_after if `sample'==1, 
				  vce(cluster st_fips)  
                                  cohort(ev_year`never') 
				  control_cohort(`c'_year)
				  absorb(tip_id##redesigntime 
				         fiscalyear##redesigntime);
  #delimit cr
  
  mat beta = e(b_iw)
  mat sigma = e(V_iw)
  erepost b=beta V=sigma, rename
  estimates save `out_dir'/sa_voucher_time_`c'_`sample'.ster, replace

 }

*-------------------------------------------------------------------------------
*--- Estimates: Event Studies All, Chain, Independent for Redesign Based on Timing of Calendar Year
*-------------------------------------------------------------------------------

cd `out_dir'

local times F5_beforetime F4_beforetime F3_beforetime F2_beforetime L0_beforetime ///
            L1_beforetime L2_beforetime L3_beforetime L4_beforetime L5_beforetime ///
	    F5_aftertime  F4_aftertime  F3_aftertime  F2_aftertime  L0_aftertime ///
            L1_aftertime  L2_aftertime  L3_aftertime  L4_aftertime  L5_aftertime 
 
foreach sample in all chain indep {
	
 preserve
 
 *Periods
 gen t=.
 local i=1
 foreach n of numlist -4/4 {
  replace t=`n' in `i'
  local ++i
 }
	
 #delimit ;
 eventstudyinteract auth `times' if `sample'==1, vce(cluster st_fips)  
                                 cohort(ev_year`never') 
				 control_cohort(`c'_year)
				 absorb(tip_id##redesigntime 
				        fiscalyear##redesigntime);
 #delimit cr	
 
 local df = e(df_r)
   
 mat b  = e(b_iw)
 mat V  = e(V_iw)
     
 *Estimates
 mat blead_before = b[1,2..4]
 mat blag_before  = b[1,5..9]
 mat beta_before = blead_before,0,blag_before
 mat beta_before = beta_before'
 svmat beta_before, names(beta_before)
 rename beta_before1 beta_before
 
 mat blead_after  = b[1,12..14]
 mat blag_after   = b[1,15..19]
 mat beta_after = blead_after,0,blag_after
 mat beta_after = beta_after'
 svmat beta_after, names(beta_after)
 rename beta_after1 beta_after
  
 mata st_matrix("se",sqrt(diagonal(st_matrix("e(V_iw)"))))
 mat se = se'
  
 mat slead_before = se[1,2..4]
 mat slag_before  = se[1,5..9]
 mat se_before = slead_before,0,slag_before
 mat se_before = se_before'
 svmat se_before, names(se_before)
 rename se_before1 se_before
 
 mat slead_after  = se[1,12..14]
 mat slag_after   = se[1,15..19]
 mat se_after = slead_after,0,slag_after
 mat se_after = se_after'
 svmat se_after, names(se_after)
 rename se_after1 se_after
 
 *95% CI
 gen lb_before = beta_before-invttail(`df',.025)*se_before 
 gen ub_before = beta_before+invttail(`df',.025)*se_before 
 
 gen lb_after = beta_after-invttail(`df',.025)*se_after  
 gen ub_after = beta_after+invttail(`df',.025)*se_after 
 
 ereturn post b V
 
 *Test 
 lincom (L0_before+L1_before+L2_before+L3_before+L4_before)/5 - ///
        (L0_after+L1_after+L2_after+L3_after+L4_after)/5

 local pv = round(r(p), 0.001)
	   
 foreach r in before after {
 
  gen att_`r'=.
  gen latt_`r'=.
  gen uatt_`r'=.

 *ATT pre period
 lincom (F4_`r' + F3_`r' + F2_`r')/3
 replace att_`r'=r(estimate) if t<-1 & t!=.
 replace latt_`r'=r(lb) if t<-1 & t!=.
 replace uatt_`r'=r(ub) if t<-1 & t!=.

 *ATT post period
 lincom (L0_`r' + L1_`r' + L2_`r' + L3_`r' + L4_`r')/5
 replace att_`r'=r(estimate) if t>=0 & t!=.
 replace latt_`r'=r(lb) if t>=0 & t!=.
 replace uatt_`r'=r(ub) if t>=0 & t!=.
 
 }

*-------------------------------------------------------------------------------
*--- Export Results: Figure
*-------------------------------------------------------------------------------

 rename t t_before
 gen t_after=t_before+.15
 
 
 if ("`sample'"=="all")   local color1 black
 if ("`sample'"=="chain") local color1 emerald
 if ("`sample'"=="indep") local color1 orange
 
 if ("`sample'"=="all")   local color2 gs5
 if ("`sample'"=="chain") local color2 dkgreen
 if ("`sample'"=="indep") local color2 dkorange
  
 if ("`sample'"=="all")   local msym1 O
 if ("`sample'"=="chain") local msym1 O
 if ("`sample'"=="indep") local msym1 T
 
 if ("`sample'"=="all")   local msym2 Oh
 if ("`sample'"=="chain") local msym2 Oh
 if ("`sample'"=="indep") local msym2 Th
 
 sum att_before if t_after>=0
 local earlyatt = round(r(mean), 0.001)
 sum att_after if t_after>=0
 local lateatt = round(r(mean), 0.001)
 
 *ssc install coefplot, replace
 *ssc install blindschemes, replace 
 set scheme plotplainblind
 graph set window fontface "Times New Roman"

 #delimit ; 
 twoway rarea uatt_before latt_before t_before if t_before<0, 
	color(gs7%50)  ||
        rarea uatt_before latt_before t_before if t_before>=0, 
	color(gs7%50) || 
	scatter beta_before t_before, mcolor(`color1') msymbol(`msym1') 
        msize(medlarge) || 
	rcap ub_before lb_before t_before, msize(vtiny) lcolor(`color1') ||
	line att_before t_before if t_before<0, lcolor(`color1') 
	lpattern(solid) ||
	line att_before t_before if t_before>=0, lcolor(`color1') 
	lpattern(solid)	||
	rarea uatt_after latt_after t_after if t_after<0, 
	color(gs12%50)  ||
        rarea uatt_after latt_after t_after if t_after>=0, 
	color(gs12%50) || 
	scatter beta_after t_after, mcolor(`color2') msymbol(`msym2') 
        msize(medlarge) || 
	rcap ub_after lb_after t_after, msize(vtiny) lcolor(`color2') ||
	line att_after t_after if t_after<0, lcolor(`color2') 
	lpattern(solid) ||
	line att_after t_after if t_after>=0, lcolor(`color2') 
	lpattern(solid) 
	yline(0, lcolor(black) lpattern(dash)) 
	xline(0, lcolor(red) lpattern(dot) lwidth(medthick))
	ytitle("ATT estimate") xtitle("Event Year")
	note("Post-EBT ATT estimate for years before package change is  `earlyatt' and for years after package" "change is `lateatt'." "P-value for difference in their average effects is `pv'")
	legend(order(3 "Before Package Change" 
	             9 "After Package Change") rows(1) position(6));
 #delimit cr

 graph export "`graph_dir'/voucher_time_sa_tip_`c'_`sample'.png", replace

 restore
  
}
}


log close
